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Abstract 

This paper is one of a series of papers describing 
the development of a new numerical method for the 
Navier-Stokes equations. Unlike conventional numeri- 
cal methods, the current method concentrates on the 
discrete simulation of both the integral and differen- 
tial forms of the Navier-Stokes equations. Conserva- 
tion of mass, momentum, and energy in space-time is 
explicitly provided for through a rigorous enforcement 
of both the integral and differential forms of the gov- 
erning conservation laws. Using local polynomial ex- 
pansions to represent the discrete primitive variables 
on each cell, fluxes at cell interfaces are evaluated and 
balanced using exact functional expressions. No inter- 
polation or flux limiters are required. Because of the 
generality of the current method, it applies equally 
to the steady and unsteady Navier-Stokes equations. 
In this paper, we generalize and extend the authors’ 
2-D, steady-state implicit scheme. A general closure 
methodology is presented so that all terms up through 
a given order in the local expansions may be retained. 
The scheme is also extended to nonorthogonal Carte- 
sian grids. Numerous flow fields are computed and 
results are compared with known solutions. The high 
accuracy of the scheme is demonstrated through its 
ability to accurately resolve developing boundary lay- 
ers on coarse grids. Finally, we discuss applications 
of the current method to the unsteady Navier-Stokes 
equations. 
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L Introduction 

This paper is concerned with the continued devel- 
opment of a new numerical method} for solving the 
Navier-Stokes equations. 1 " 3 The distinguishing fea- 
tures of the current method have been previously de- 
scribed, and may be summarized as follows: The 

current method (i) provides for a unified treatment 
of space and time; (ii) represents the local discrete 
primitive variables through Taylor series expansions 
that identically satisfy both the integral and differ- 
ential forms of the Navier-Stokes equations; (iii) bal- 
ances fluxes across cell interfaces as an integral part 
of the numerical formulation; (iv) evaluates fluxes at 
cell boundaries using exact functional expressions (to 
the order of accuracy of the local expansions); and (v) 
solves explicitly for the unknown derivatives of the lo- 
cal discrete primitive variables. 

Scott and Chang 1 developed a Newton’s method 
implicit scheme for the two-dimensional, steady 
Navier-Stokes equations. Calculations of laminar flow 
in a channel showed that the developing boundary 
layer could be accurately resolved using as few as six 
cells per channel width, Scott 3 extended the internal 
flow scheme to external flows, and showed that the 
boundary layer on a finite flat plate could be accu- 
rately resolved using only ten cells from the wall to 
the free stream. Chang 2 developed an explicit scheme 
for the one-dimensional, unsteady Navier-Stokes equa- 
tions, and presented numerical results for the shock 
tube problem with viscosity. It was shown that the 
shock and contact discontinuities could be crisply re- 
solved without the use of flux limiters or weighting 
functions. 

In this paper we are concerned with the fur- 
ther development of the steady-state flux conserva- 
tion scheme first presented in Reference [1]. Here we 

}Known most completely as The Space-Time Conser- 
vation and Solution Element Method, and frequently 
abbreviated as the STS Method, Conservation Ele- 
ment Method, Solution Element Method, or the CE- 
SE Method. 
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generalize the scheme by (i) providing a general clo- 
sure methodology so that all discrete unknowns of the 
same order in the local expansions may always be re- 
tained and (ii) by extending the scheme to nonorthog- 
onal Cartesian grids. Calculations of a variety of flow 
fields indicate that developing laminar boundary lay- 
ers can be accurately resolved on grids which are much 
coarser than those employed by conventional numeri- 
cal methods. 

In the next section we begin by considering the 
general problem of steady-state flux conservation in 
two space dimensions. Following that, we summarize 
the steady-state flux conservation scheme as originally 
presented in [1] and [3]. Then we show how the scheme 
may be generalized and extended to general flow fields 
and nonorthogonal grids. Finally, we discuss applicar 
tions of the present approach to the unsteady Navier- 
Stokes equations. 


solution to boundary value problem (2.2). In addition, 
let eaci hm be analytic throughout V . (From here on, 
we drop the m subscript for clarity.) Then near any 

(x 0 ,y 0 ) in V, hm h — (h r , h y ) may be expressed by 
the uniformly convergent Taylor series expansions 3 

h*{x,y) = h*(x 0 ,y 0 ) + ^-(i 0 , y 0 ) (* - * 0 ) 


+ + - (2.3) 




(r - *„)"'* (y - y,)* 


n=0*=0 


*dy k 


(n - k)\ 


k) 


h v (x,y) = h*(x 0 ,y 0 ) + 


dh* 

dx 


( x 0 ,y 0 )(x-x 0 ) 


II. Steady-State Flux Conservation 
in Two Spprg pi mensions 

We begin by considering the system of integral 
conservation laws 


/ h m Ts = 0 (2.1) 

JsiV) 

in Ez, where h m , m = is a steady-state flux 

density vector which is a function of M intermediate 
(p rimit ive) variables. S(V) is an orient able surface 
inclosing an arbitrary volume V in two-dimensional 
space, and ds is the outward area vector which is nor- 
mal to S'(V’). Equation (2.1) is the governing equation 
for steady-state flux conservation in two space dimen- 
sions. 

Suppose that (2.1) holds for any V in a fixed 
domain V in £3, and that each h m varies smoothly 
throughout V. Then one obtains, by way of the diver- 
gence theorem, the differential conservation law 

Vh m = 0, (2.2) 

where V = (^, 5^). If information about the be- 
havior of the h m (or their intermediate variables) on 
the boundary B(V) of V is also known, then equa- 
tion (2.2), along with the known boundary informs 
tion, represents a steady-state boundary value prob- 
lem. We will subsequently refer to equation (2.2) and 
its associated boundary conditions as boundary value 
problem (2.2). 

Suppose that V is the rectangular region shown 
in Figure 1, and that h m , m = is the exact 


dh v ( w x 
+ ^-(«o,yo)(y-yo) -1- ... 


00 n 

= ££[ 


n=0*=0 


dx n ~ k dy k 


(*o,yo)] 


(2.4) 

(x - x 0 ) n-fc (y - y 0 ) k 
(n - Jfc)! k\ ’ 


If the infinite series (2.3) - (2.4) is truncated at some 
point, then it becomes an local approximation to the 
exact solution h. 

Let V be discretized by the rectangular mesh 
shown in Figure 2, and let (a? s - , yj ) denote the center 
point of each cell. Let (£2)1,; denote the order-two 
Taylor series expansion of h about the point (x;, t/j). 
Then 

(*5)W = (2-5) 


V'V' 'f 9 n h z / \1 (x — x,) n * 

- J (n — jfc)! 


(y-yj)* 

Jfc! 


and 


= (2-6) 

- ffi (x v .<\ (izfiUlt 

- 2-,2-s [dx n - k dy k( " Vl) l (n — jfc)! Jfc! ' 


n=OJk=0 


As an approximation to h on the (*‘,j)th cell, it 
can be shown that 3 the error of (hjJij is bounded 
by ^*[max(Ax, Ay)] 3 , where Mij is a local bound 
on the third derivatives of h. Thus, {(£2)14 = 

1,..., Niyj = 1, ..., JSTj } is a discrete approximation to 
the exact solution h on X>. 

The local Taylor series expansion (£2)14 is a so- 
lution of the differential conservation law (2.2), since 
its divergence is identically zero. 3 This follows from 
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the fact that (2.3) - (2.4) also satisfies equation (2.2). 
Since (h2)ij is a second-order polynomial that pro- 
vides the exact values of h and its derivatives through 
order two at the point (a?,*, jy), it will be referred to as 
the exact second-order solution of equation (2.2). 

Since the derivatives of h are well-defined 
throughout V, we may now consider the approximate 
second-order solution h 2 to boundary value problem 
(2.2), where 


with 


h 2 ^ h = ( h*,h y ) (2.7) 

h*(x,y;i,j) *= (2.8) 


Ko + *"•(*“**) + A «.i (y-w) + 
Ko( x ~ x i ) 2 + *t»(* “**)(»-») + KAy-ytf 

and 


h v (x,y,i,j) ^ (2.9) 

A*. + A* 0 (x — Xi) + A» tl (y -«,•) + 

ftf >0 (x — *j) 2 + A» 1 (ar — *<)(y — «#) + A* a (y - y,) 2 . 

The coefficients A* 0 , **„, etc., are the numerical ana- 
logues of h*(xi,yj), ^r(xi,yj), etc., and are un- 
knowns that are to be determined (the i ) j subscripts 
have been omitted for clarity). Generally, the un- 
known coefficients above are expressed in terms of one 
or more intermediate variables. This will become clear 
later, but need not concern us here. 

To faithfully represent h throughout V> the dis- 
crete approximation {(h)ij , * = h ^**3 = •••» fy} 

should: (i) satisfy (2.1) on each cell and each union of 
cells in V\ (ii) have zero divergence; and (iii) satisfy 
the boundary conditions on B(V). 

Let the vertices of the (i, j ) th cell be denoted by 
and 5 as shown in Figure 3. Then (i) above 
requires that 

<f hTs = 0, (2.10) 

JPQRSij 

where we have taken the integration to be positive in 
the counterclockwise sense. Let [J(PQ)]ij denote the 
flux of h through the line segment PQij , and similarly 

for J{QR), J(RS)i and J{SP ). Then, omitting the 
i y j subscripts, one obtains the flux expressions 

J(PQ) = (2.H) 


Ax 3 

12 


h y 2Q + Ax + A S.o] 


Ay 3 ,, 

12 h '- 3 


J(QR) = 

Ax 2 
4 


( 2 . 12 ) 


Ay A* 0 


A^ 1 * I LX 

2 A x ,8 ' *0,0 


J(RS) = 


Ax 3 

12 


A?,. - 


[ Ay 2 


A* PJ-A?., 


(2.13) 

A y^y + /,» 1 

2 '*0.1 ' O.oj 


+ A, [^-M, + %K. + <•]• 

Equation (2.10) requires that the fluxes satisfy 
J(PQ) + 7(Ql) + J(E3) + J(SF) = 0. (2.15) 
Thus, one obtains the flux conservation constraint 

h*, 0 + Kr = 0. (2-16) 

In view of (2.10), each cell in the mesh is called an 
conservation element. Since h may be discontinuous 
across cell interfaces, each cell is also called an solution 
element . 4 

Imposing condition (2.16), one obtains from 
(2.11) - (2.14) the following normalized flux expres- 
sions: 

J ( P Q1 = (2.17) 

Ax 

Ax 2 , y Ay 2 , y Ay , c , y 

^2"A* 0 + —K a - — A 1(0 + *»,„ 


J{SP) = 
rAx 2 


(2.14) 


J(QR) 

Ay 


(2.18) 


Ay 2 , * Ax 2 x Ax . x 

^2 *o,a + ^ *a,o 2 ftl, ° ‘ n °’° 


/(JJS) 

Ax 

Ax 2 , # Ay 2 , y Ay , x , y 

-j2~ A?. 0 + — A;,a + -y A 1>0 + A o,o 


(2.19) 
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J(5p) 

Ay 


( 2 . 20 ) 


Ay 2 

12 


h 


X 

0,3 


+ 




+ 


Equations (2.16) - (2.20) are sufficient to ensure 
that (2.1) is satisfied on each cell. To ensure that (2.1) 
is satisfied globally, it is necessary that fluxes be con- 
served across cell interfaces. 

Let the (i,j) th conservation element be denoted 
by CE(i,j). Consider the vertical interface between 
CE(iJ) and CE(i + l,j). The flux leaving CE(iJ) 
through this interface is given by [/(5P)], j. On the 
other hand, the flux leaving CE(i + 1 , j) through this 
same interface is given by [J((j5)]i+ij. Flux conser- 
vation across the interface requires the sum of the two 
to be zero. Thus, one obtains the interface condition 


(2.24) ensure the satisfaction of a necessary condition 
for the second-order discrete unknowns to converge to 
their respective analytical counterparts. 

We should also emphasize that (2.23) - (2.24) are 
local conditions, as is the first-order constraint (2.16). 
Thus, a high-order-accurate simulation of the differ- 
ential conservation law is possible without stringing 
together more and more mesh points. 

The above approach provides the conceptual 
framework for the discrete simulation of steady-state 
flux conservation in two space dimensions. We now 
demonstrate the applicability of the framework to the 
Navier-Stokes equations. 

HI. Conservation of Mass and Momentum 
in Incompressible, Viscous Flow 

A. Discrete Flux Vectors 


Ay 2 ,, Ax 2 Ax 

^2 ^°» a 4 ^2,o ' 2 ***♦• ' ^°*° 




Ay 2 ^ , Ax 2 


Ax , 


22 4 ^s,o 2 ^ 1 *° 


( 2 . 21 ) 


= 0. 


J »+i j 

Similarly, across horizontal interfaces, one obtains the 
interface condition 


h* + - —h* + h 9 

22 f * a «° ' 4 '*0,3 2 Xt ° T 0,0 


( 2 . 22 ) 


j tj 


= 0 . 

J y+i 

Equations (2.16) - (2.22) above ensure that the 
integral form of the governing conservation law will be 
satisfied throughout V. 

Having provided for the satisfaction of the more 
fundamental integral conservation law, we may now 
turn our attention to the differential conservation law 
(2.2). Requiring h to have zero divergence on each 
conservation element, one obtains the two second- 
order constraints 


—v> 
12 ao 


+ **• + **., 


2 A*. + h 9 tl = 0 (2.23) 

2 hi, + = 0. (2.24) 

By virtue of equations (2.23) - (2.24), the approximate 
second-order solution h satisfies (2.2) identically, as 
does the exact second-order solution &2* In Reference 
[3] it was rigorously proved that conditions (2.23) - 


The governing equations for the conservation of 
mass and momentum in a steady, two-dimensional, in- 
compressible flow field can be written 5 



SIS’ 

+ 

<§>!§> 

ii 

o 

(3.1) 

d , 2 

. d , , 

+ P - T„) + ^(ut> - T *y) - 0 

(3.2) 

d . 

Tx (uv 

~ Try) + ^(V 2 + P - Tyy) = 0 

(3.3) 

where 

r _ 2 du dv 

xx dx dy 

(3.4) 


_ 1 ( du dv . 

Txy Re L + dx 

(3.5) 


T — —2— (2 — _ 

Tyy ~ SRe^dy dx ] ' 

(3.6) 


The independent variables x and y are horizontal and 
vertical Cartesian coordinates, respectively, and u and 
v denote their respective velocity components, p de- 
notes the static pressure, and Rcl is the Reynolds 
number, -, where the viscosity p and density p 
are assumed to be constant. L and U 0 0 refer to some 
referei ce length and velocity, respectively. 

The integral form of (3.1) - (3.3) is given by 


® h„ ■ ds = 0 

(3.7) 

/S(V) 

P h xu ds - 0 

(3.8) 

S(V) 
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( 3 . 21 ) 


{ h YU ds = 0 (3.9) 

Js(V) 

where h u , h X u> and h YM are defined by 

h„ =' (u,v) (3.10) 

h X M — (t* 2 + P — i uv — r *y) (3.11) 

h YK ( d — (tit) — Try, V 2 + p — Tyy). (3.12) 

Equations (3.7) - (3.9) are a coupled system of in- 
tegral conservation laws whereby mass and momentum 
are conserved by way of the primitive (intermediate) 
variables tx, v, and p. The system of equations (3.7) - 
(3.9) is of the same form as system (2.1). 

Corresponding to the continuous flux vectors de- 
fined by (3.10) - (3.12^, we introduce the second-order 
discrete flux vectors and h Yhi , where 



d - («>«) 

(3.13) 

- dej 

Qxm — 

(tX 2 + p — T xx , tit) - Z T v) 

(3.14) 

h & 

*}ym — 

(u V T X y , V -f* P Tyy ) 

(3.15) 

and 



Tr* = 

^-(2du/dx - dv/dy) 

(3.16) 

T xy = 

= — (dv/dy + dv/dx) 

(3.17) 

ii 

tJ-(2 dv/dy - dv/dx). 

oltCi, 

(3.18) 


By definition, the discrete primitive variables tx, v, and 
p are then second-order expansions which may be 
written 

v(x,y;i,j) =/ (3.19) 

tio.o + «!,,(*-*<) + «o.i(y-y;) + 

u a>0 (* - xi) 2 + -*<)(»-») + «».*(y-%) 2 

v(x,y;i,j) d - (3.20) 

v 0 ,o + t)i, 0 (x-Xi) + t) 0tl (y-y,) + 
v a , 0 (x-xi ) 2 + — x,)(y—yj) + v 0 , a (y-yj ) 2 


P(x,y;i,j) 

Po.o "b Pi,o( x — *•) d" Po.i(y SO') “b 

P2,o( x ~ x i) 2 + PiA x - x i)(y - vs) + PoAv-yjf- 

For clarity, the i,j subscripts have been omitted from 
the discrete Taylor coefficients ix 0t0 , v 0t0 , p 0 ,o> etc. 
These coefficients are the unknowns to be solved for. 

Each of the discrete flux vectors h XMy and 
h YM may be expressed in the form of (2.7) - (2.9), 
where h* 0 , Aj |0 , etc., are functions of tx 0>0 , v 00 , p 0(0 , 
tx 1>0 , t> 1>0 , p lf o, etc. For example, in reference to h = 
h XM} h* 0 is the constant term of the expression u 2 -h 
p - r xx , which is given by tXo,o+Po,o- 3 ^ 7 ( 2 ui,o*-t;o.i). 
In the Appendix, we present the functional expressions 
for each component of the discrete flux vectors h Ml 
h X jXi aud 

B. Constant-Pressure-Gradient Flows 

1, Laminar Flow in a Straight Channel 

Consider the channel geometry and associated 
mesh shown in Figure 4. On each solution element, 
we represent tx, v, and p through the second-order ex- 
pansions (3.19) - (3.21). There are 18 discrete un- 
knowns on each cell. Since there are 18iVj Nj unknowns 
altogether, 18 NiNj conditions are required to have 
a closed system of equations. Thus, for the present 
equilibrium-type boundary value problem, closing the 
system requires a global matching of numbers of con- 
ditions and unknowns. 

The conditions to be imposed are obtained from 
Section II. For each of the three conservation laws 
(3.1) - (3.3), we must impose the local constraints 
(2.16), (2.23) and (2.24) on each conservation element. 
This provides 9 NiNj conditions. (See equations (A.l) 
- (A.9) in the Appendix.) The interface conditions 
(2.21) and (2.22) must also be satisfied. This pro- 
vides an additional 3 Nj(Ni — 1) 4* 3 Ni(Nj — 1) = 
6 NiNj - 3 Ni - 3 Nj conditions. 

Another iNi+ZNj conditions are obtained by way 
of boundary conditions. At the channel inlet, we spec- 
ify the two components of velocity, and at the exit we 
specify the pressure. Along the upper and lower walls 
we require that there be zero mass flux. In addition, 
we impose the no slip condition for tx at the midpoint 
of the wall face of each cell. 

For laminar flow in a straight channel, it is well 
known that the streamwise pressure gradient is nearly 
constant, except in a short entrance region. Thus, 
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« 0 and « 0 throughout most of the channel. 
It is reasonable then to make the simplifying assump- 
tion that p 3t0 = 0 and p ltl = 0 on each cell. This 
provides an additional 2N{Nj conditions. 

The number of conditions still needed to close the 
system is then Ni(Nj — 1). Since this is the num- 
ber of horizontal interfaces in the mesh, there is an 
additional degree of freedom in specifying horizontal 
interface conditions. We choose to close the system 
by requiring the u velocity to be continuous at the 
midpoint of each horizontal interface. 

The discrete boundary value problem outlined 
above is a coupled system of second-order polynomial 
equations in the unknowns u 0p0 , v 0tQ , p 0t 0 , etc. How- 
ever, six of the higher-order unknowns can be eas- 
ily eliminated from the system using the local con- 
straints (A.l) - (A.9). (See equations (A.10) - (A.15).) 
The total number of unknowns that must be explicitly 
solved for is thus reduced from 16N{Nj to 10 NiNj. 
Of the equations that remain, 3 NiNj represent local 
conditions, and the remaining 7 NiNj represent inter- 
face and boundary conditions. Numerical solutions to 
the nonlinear system of equations may be obtained by 
Newton’s method. 1 The Jacobian matrix has the form 
shown in Figure 5. 

In Figures 6 and 7 we present numerical re- 
sults from calculations performed at Reynolds num- 
bers (based on channel height h) of 100 and 2000. 
The specified inlet velocity profile for each case is 
u\ = | [l — [2 (y— |)] 2 * * * 6 * ] (following Dill and Himansu 8 * ). 
Along with the numerical results, we also show the in- 
let profile and fully developed analytical solution. The 
Re = 100 results were obtained from the 39 x 6 grid 
shown in Figure 4. Due to the thinner boundary layer, 
the Re = 2000 results were obtained from a 39 x 10 
mesh. For each case, Newton’s method converged to a 
maximum residual error of 10" 10 in four and five iter- 
ations, respectively (starting from uniform flow). The 
corresponding CPU times on a Cray YMP were 1.1 
and 3.5 seconds. 

2. High- Reynolds- Number Flow 

Around a Finite Flat Plate 

In Figure 8 we show a typical mesh for calculation 

of the thin airfoil boundary layer problem. The airfoil 

lies on the x-axis between 0 and 1, and the mesh in- 

cludes both an upstream and wake region. The grid 

spacing is exponentially stretched in the y direction, 
and the spacing may also be nonuniform in the x di- 
rection. 

As in the channel flow problem, there are a total of 

18 NiNj unknowns to be solved for. The methodology 


for closing the system here closely follows that given 
in the previous section. 

Let N a denote the number of solution elements 
between the airfoil leading and trailing edges. Then 
(2.16) and (2.21) - (2.24) provide 9 NiNj + 3 N s (Ni - 
l) + SN i (N j -l)-ZN a = IbNiNj-ZNi-SNj-SNa 
conditions. 

Boundary conditions provide an additional 4 Ni + 
3 Nj + 4 N a conditions. For each cell adjacent to the 
airfoil, we require the mass flux through the wall face, 
and the ti velocity at the midpoint of the wall face, to 
be zero. At the upstream boundary, we specify the ve- 
locity, and downstream we specify the pressure. Along 
the free-stream boundary cells, we specify zero y gradi- 
ent conditions for ti and v. Imposing the zero pressure 
gradient condition, we set p 2 , 0 = 0 and p lpl = 0 on 
each cell, for an additional 2 NiNj condtions. Finally, 
an additional Ni(Nj — 1) — N a conditions are obtained 
by requiring the ti velocity to be continuous at the 
midpoint of horizontal cell interfaces. 

In Figure 9 we compare the predicted boundary 
layer profile with the Blasius solution at four differ- 
ent locations along the airfoil. The Reynolds number 
based on airfoil chord is 100,000. The results were ob- 
tained from an 81 x 20 grid with uniform x spacing and 
nonuniform y spacing with 12.5% exponential stretch- 
ing. Only ten cells are used across the boundary layer, 
since the flow field is computed on each side of the air- 
foil. The y spacing at the wall was A y w = .0015, and 
the free stream boundaries were located at y = ±.027. 
The CPU time required on a Cray YMP was 50.7 sec- 
onds (nine Newton iterations), with uniform flow as 
the starting solution. 

Figure 10 compares the predicted v velocity at x 
= .805 with the Blasius v velocity. The differences 
that appear near the edge of the boundary layer are 
believed to be due to finite airfoil effects [7, p. 137]. 

Ie Figure 11 we show the pressure coefficient in 
the trailing edge region. The results were computed on 
an 110 x 22 grid that was refined near the trailing edge 
with exponential x stretching. The x spacing at the 
trailing edge was .007. Our results show the pressure 
coefficient to have a minimum value of - .014, which 
agrees well with the Reduced Navier-Stokes (RNS) cal- 
culations of Srinivasan and Rubin. 6 

C. Ge neral Flow Fields 

1. Generalized Closure Procedure 

For a general two-dimensional problem with an Ni 
x Nj mesh, the local expansions (3.19) - (3.21) provide 
a total of IS NiNj unknowns. Thus, IS NiNj condi- 
tions are always required to close the discrete system of 
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equations. In Section II. B., we imposed UMVfJVj gen- 
eral conditions, and then added the 2 N{Nj special con- 
ditions \p 2$0 ]i j = 0, \pi,i]ij = 0. Since the general con- 
ditions already ensure that the integral and differential 
forms of the governing conservation laws are fully sat- 
isfied , it will be necessary to introduce higher-order 
conditions to close the system. Such higher-order con- 
ditions can be imposed either on the discrete primitive 
variables, or on the discrete flux vectors. Due to the 
more fundamental nature of the flux vectors, the sec- 
ond alternative is generally preferred over the first. 

We are thus led to consider the first-order dipole 
moment interface conditions 


and 


f (x - Xi)(hij - hij+i) • ds = 0 

Jsx 

f (y- - Ei+ij) ds = o, 


(3.22) 


(3.23) 


where x\ and x 2 denote the endpoints of an horizon- 
tal interface between C E(i t j) and CE(i,j + 1), and 
yi and denote the endpoints of an vertical inter- 
face between CE(i, j) and CE(i+lJ). Since aerody- 
namic flow fields are generally dominated by gradients 
in the cross-flow direction, equation (3.22) will usually 
be preferred over equation (3.23). However, both con- 
ditions are fully justifiable, since the exact solution h 
satisfies both equations. 

For a rectangular grid, equations (3.22) and (3.23) 
become 


K+f **.]„ - - 0 < 3 - 24 > 


Ay, 


and 


<.]„ - (325) 


respectively. 

We now consider again the channel and flat plate 
flow fields. For channel flow, we impose (3.24) for 
h = h X M and h = h YM , thereby obtaining an addi- 
tional 2 Ni(Nj — 1) conditions. The final 2 Ni conditions 
may be obtained by imposing am additional boundary 
condition on the upper and lower wadis. From the y- 
momentum equation one obtadns the pressure bound- 
ary condition 5 

dp _ 1 d 2 v 

dy ~ Itedy 2 ' 

which may be imposed at the center point of the wall 
face of each cell. (Equation (3.26) is not a boundary 



condition in the usuad sense, but is instead a differ- 
ential equation to be satisfied on a wall parallel to 
the x-axis. Even though we already require the y- 
momentum equation to be satisfied identic adly on each 
cell, equation (3.26) still represents an independent 
condition which may be imposed adong the wall.) 

For the flat plate problem, imposing (3.24) for the 
x- and y-momentum equations provides 2 Ni(Nj — 1) — 
2 N a conditions, where N a is the number of cells be- 
tween the leading and trailing edges. Applying (3.26) 
to each side of the airfoil gives 2 N a conditions, and im- 
posing the zero pressure gradient condition as a free- 
stream boundary condition provides the final 2 N{ con- 
ditions. 

In general, conditions (3.22) - (3.23), together 
with additional boundary conditions, cam be used to 
close the discrete system of equations for an arbitrary 
flow field. By introducing such higher-order condi- 
tions, the system is not only closed, but additional 
physics are brought into the scheme. Thus, it is pos- 
sible to adapt one’s particular choice of closure to the 
physics of a specific flow field. 

In the general case, there are an additional 2 NiNj 
unknowns to be solved for. However, the size of the 
system to be solved is still 10 NiNj, as in the con- 
stamt-pressure- gradient scheme. The reason is that, 
in the general case, the two new unknowns p 3 , 0 and 
p lx cam be explicitly eliminated from the system. (See 
equations (A.20) - (A.21).) Thus, a total of SNiNj 
unknowns are eliminated, and 10 NiNj still remain. 

Application of the generalized scheme described 
above to the channel flow problems of Figures 6 and 
7 reproduces the earlier results. The comparison is 
shown in Figures 12 - 13. 

For flat plate flow, the neglected terms p 3 , 0 and 
Pi a become significant for Reynolds numbers less 
than about 7500. Figures 14 and 15 compare re- 
sults from the two schemes at Reynolds numbers 
of 7500 and 1000. As Re decreases, the boundary 
layer becomes much thicker and the pressure gradi- 
ent varies significantly throughout the flow field. At 
the lower Reynolds number, the improved accuracy of 
the general scheme over the constamt-pressure-gradient 
scheme is readily apparent. 


2. Extension to Nonorthogonal 
sjfm Grids 

Extension of the generalized scheme presented 
above to nonorthogonad grids can be accomplished in a 
straightforward manner. Consider the nonorthogonal 
mesh shown in Figure 16. We denote the mesh points 
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by = 0 and the cen- 
ter points by = 1. -i NiJ = 1. •••> N j}> 

where x ei j — (z.'-ij-i + *tj-i + *ij + *»-ij)/ 4 

and y eiJ + Vij-1 + Vij + lfc-ij)/ 4 - 0“ 

each solution element, the discrete flux vectors and 
primitive variables are expanded about the cell center 

Since we require h Uy hxMi and h Y ht to be identical 
solutions of the differential conservation law (2.2), lo- 
cal flux conservation on each CE(i f j) is assured. This 
follows from the divergence theorem. Consequently, 
to guarantee global flux conservation, we need only 
ensure that fluxes are conserved across the nonorthog- 
onal cell interfaces. This in turn requires that the nor- 
malized flux expressions (2.17) - (2.20) be generalized 
to account for the nonorthogonality of the mesh. 

The generalized fluxes are derived most easily us- 
ing a local coordinate system relative to the cell center 
(see Figure 17). The cell vertices P y Q , JZ, 5, relative to 
(scij.yaj), are given by (^p,^p), 

(^,-4*-), where 


A X p def 




(3.27a) 


4- A 1 (Ci f o + *nCo t i 4- SCi t \ + 2m 5 Co, 2 ) 
+ A 2 (C2 f o + ™C\ t i + m 2 Cqj) 


where 


m 'i' (^ - ^ ) (3.29a) 

6 'W ^ l + (3.29b) 

a 1 tg i ( d*, _ a^, (3Mc) 


J(RS) 

Ax T t At, 
2*2 


a 1 ‘a + (3.29d) 


= Co,o 4- £ Co,! + 6 2 Co t 2 (3.30) 


4- A 1 (Ci ? o 4- m Cq,i 4- 6 Ci t \ + 2 mS Co , 2 ) 
+ A 2 (C*2,o + m Ci,i + m 2 Co, 2) 


Ay p del 
— 2 - 

(3.27b) 

Ax. def 
2 ~ Xci J 

(3.27c) 

Ay, def 
2 — V*-ij — 

(3.27d) 

Ax r del 
2 — Xci J “ 

(3.27e) 

Ayr def 

~2~ — Vcij—tti-ij-x 

(3.27f) 

Ax, def 

~2~~ — x m’- 1 ” x «**i 

(3.27g) 

Ay, de/ 

~ 2~ ~ ycij — 

(3.27h) 

We assume that etc., are 

all positive. 

For the flux through an PQ or RS interface, let 

C„-jfc,t — Ki-h,k ~ 

for n = 0, 1, 2, fc = 0, ..., n, 
where m will be defined shortly. Then 

(3.28) 


where 


def /A y r Ay, x ./Ax r Ax,> 

m - <— - 2 )/( 2 + 2 > 

(3.30a) 

def Ay r Ax r 

v — — ^ 4* m _ 

2 2 

(3.30b) 

* 1 ^ 5 / ^ / Ax, Ax r ■, 

2^2 2 

(3.30c) 

and 


*2 1 [( Arr ) 2 A*r Ax, ( ^Ax,j 2 

3 L 2 2 2 2 

] (3.30d) 

The normalized first-order dipole moment flux 
through an PQ interface is 

11 

ni 

(3.31) 


A 1 (Co,o + £Co,i 4- $ 2 Coj) 

4- A 2 (Ci,o 4* mCo,i 4- £Ci f i 4- 2 m 6 60 , 2 ) 


j (PQ } ~ = C 0 ,o + 6C 0 ,t + 6 2 Co, 2 (3.29) 


4“ A 3 (C 2 ,o 4- wiCij 4- m 2 Co, 2 ) 


where 
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A 3 d = (3.31a) 

1 r / Axp \ 3 , Ax P s i Ax q Ax p( Ax t ,2 /Ax^ai 

4 J — + —V— ; V 2 > J* 

and all other parameters are defined as in (3.29a) - 
(3.29d). Similarly, the normalized first-order dipole 
moment flux through an RS interface is 


/(AS) 

Ar r , Af 
„ 2 “ r 2 



(3.32) 


A 1 (Co,o + SCo t i + ^ Co,a) 

4 A 2 (Ci,o + mCo,i 4* £Ci f i + 2m6Cv t 2) 

4 A 3 (C*2 } o 4" rn Ci,i 4" m 2 Co, 2 ) 

where 

A 3 (3.32a) 

1 r/Aar #x 3 / Ax, x 2 Ax r ( Ax, f Az r ^2 /Ax rN 3l 

4 ^~ 2 ~’ ~ 2 ~ + ~ 2 ~ [ ~ 2 ~ ) 2 ' J’ 

and all other parameters are defined as in (3.30a) - 
(3.30d). 

For the flux through an QR or SP interface, let 
C„_ M = K_ ktk - mh» n _ ktk (3.33) 


for n = 0, 1, 2, k = 0, n. 
Then similar to the above, we have 


— a — = ^ 0,0 + *Cw 4* ^ 2 C*2,o (3.34) 

4- A 1 (Co,i 4 mCi f o 4 SCi t i 4 2m6C t 2,o) 

4 A 2 (Co,2 4 roCi,i 4 ro 2 C2,o) 

where 


m 


dej AXr 

2 


Ayr 

2 


a 1 <(£/ i 

2 V 2 2 ’ 


and 


(3.34a) 

(3.34b) 

(3.34c) 


A 2 W | [(^) 2 - + (^) 2 ] (3.34d) 


Ay» , Ay, 


where 


m 


dej Ax, 

^ = 1“ + 


2 v 2 


and 


1 4 £ 2 C2,o 

(3.35) 

li + 2mSC2,o) 
(- m 2 Cj.o) 

.+^t) 
2 7 

(3.35a) 

Ay, 

2 

(3.35b) 

Ay,x 
2 ' 

(3.35c) 


(3.35d) 


The normalized first-order dipole moment flux 
through an djR interface is 


J(QR) 




J dp 


(3.36) 


A 1 (Co,o 4 6 C\ 1 0 4 6 2 C 2 ,o) 

4- A 2 (Co,i 4 mCi,o 4 £Cj,i 4 2m£C2,o) 

4- A 3 (Co , 2 4 niCi,i 4 m 2 C2,o) 

where 

A 3 =' (3.36a) 

1 f/Ay n 3 /Ay n 2 Ay r Ay q( Ay r ^i ,Ay r v3' 

and the remaining parameters are defined by (3.34a) - 
(3.34d). 

Finally, the normalized first-order dipole moment 
flux through an SP interface is 


J(SP) 


L 2 ^ 2 


dp 


(3.37) 


A 1 (Co,o 4 6 C\ t o 4 5 2 C 2 ,o) 

4 A 2 (Co,i 4 m Ci'O 4 $C\ t i 4 2m5C2,o) 
4- A 3 (Co , 2 4 m Ci r i 4 m 2 C 2 ,o) 
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where 


A 3 d = (3.37a) 

\ [(^) 3 - (^)’% + ^(%)’ - (^) s ] . 

and the remaining parameters are defined by (3.35a) - 
(3.35d). 

To complete the nonorthogonal formulation, one 
must modify the pressure boundary condition (3.26) 
to account for any variation in the wall slope. The re- 
maining boundary conditions require no special treat- 
ment. 

In Figures 18 - 20 we consider channel flow with 
a converging or diverging section. For each case, the 
ramp length is one channel height, and the inlet ve- 
locity is Uj = | [l — [2 (y — |)] 6 ] . The ramp angles 
for the three cases are 10, 15, and 5 degrees, respec- 
tively. In Figures 18 b. and 20 b., we consider flow at 
a Reynolds number of 100. The calculations were run 
on the 32 x 8 grids shown in Figures 18 a. and 20 a. 
The two cases converged in six and eight Newton itera- 
tions, respectively. For both cases, the results indicate 
that the flow becomes fully developed by the end of 
the channel. In Figure 19 b., results are presented for 
flow at a Reynolds number of 200. The calculations 
were performed on the 32 x 10 grid shown in Figure 19 
a. For this nonsymmetrical flow at a higher Reynolds 
number, Newton’s method required ten iterations to 
converge (11.0 seconds on a Cray YMP). Unlike the 
results in Figures 18 b. and 20 b., the flow does not be- 
come fully developed by the end of the channel. This 
is consistent with the theory for developing laminar 
flow in a channel [7, p. 186]. 

A similar calculation to the above, using a general 
coordinates STS implicit scheme, has been carried out 
by Dill and Himansu. 8 

D. Accuracy of Derivatives 

The numerical results presented thus far have 
demonstrated that the current method can resolve de- 
veloping boundary layers on coarse grids. This clearly 
indicates that the discrete derivatives, at least through 
some order, are being obtained accurately. The accu- 
racy of the discrete derivatives is a current research 
topic, and will be discussed in detail in a future paper. 

Our experience to date indicates that the accuracy 
of the first-order derivatives is in general very good. 
On the other hand, the accuracy of the second-order 
derivatives varies considerably, and appears to depend 
on a number of factors. 

As a model problem for estimating the accuracy 
of the derivatives, we consider the flow field around a 


finite flat plate. The Blasius similarity solution pro- 
vides the first-order asymptotic solution to the Navier- 
Stokes equations which is valid near the wall and away 
from the leading and trailing edges. 9,10 Thus, at very 
high Reynolds numbers, the Blasius solution can be 
used to obtain analytical values for the derivatives of 
ti and v. 

In Figure 22, we compare the predicted (i.e., 
u 0t i) with the Blasius There is in general very 
good agreement between the two. (The numerical is- 
sues involved in obtaining accurate derivatives of the 
Blasius solution are addressed in a paper which is cur- 
rently being written. 11 ) The large value of |“ near the 
wall is also evident. 

Figure 23 compares the predicted (i.e., 2ti 0 , 3 ) 
with the Blasius |^y. Although not as accurate asu 0)1 , 
the relative error of u 0|3 is of the order of 10“ 3 at a 
number of grid points. Note that the magnitude of 
is greater than 10 3 throughout mo6t of the boundary 
layer, and obtains a maximum value which is greater 
than 10 4 . 

IV. Space-Time Flux Conservation for the 
Unsteady Navier- Stokes Equations 

The discrete formulation presented above for 
the steady, incompressible Navier-Stokes equations 
is designed and constructed to simulate steady-state 
boundary value problems which are of equilibrium 
type. However, when considering the unsteady Navier- 
Stokes equations, one is confronted with an initial- 
boundary value problem which is of a marching type. 
Thus, the problem of space-time flux conservation es- 
sentially reduces to how to march the discrete solution 
forward in time in such a way that space-time flux 
conservation is always satisfied. In general, one may 
advance the solution in am explicit or implicit man- 
ner. In either case, the focus of the current method 
is to ensure that fluxes are conserved in space-time, 
through the use of Taylor series expansions that iden- 
tically satisfy the governing integral and differential 
conservation laws. 

As a model problem for the unsteady Navier- 
Stokes equations, Chang 2 considered the shock tube 
problem with viscosity. The governing equations are 
the 1-D, unsteady Navier-Stokes equations. Let p, v , 
p, and 7 be the mass density, velocity, pressure, and 
constant specific heat ratio, respectively. Let 

«i = p (4.1) 
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«2 = pV 

«3 = p/(y- i) + (i/2V>t> 2 

h =«2 

h = (7 - 1)«3 + (l/ 2)(3 - 7) («2) 2 /«i 
/ s = 7 « 2 «a/«i - ( 1 / 2 )(t - l)(« 2 ) s /(«i ) 2 
and let 

A = 0 

r 4_ «2 

/2 ~ 3i2e m 

7 r «3 _ («2) 2 1 

RePr l«i 2 («i) 2 J ’ 


h = 


3 fie \ui 


B 


(4.2) 

(4.3) 

(4.4) 

(4.5) 

(4.6) 

(4.7) 

(4.8) 

(4.9) 


where Re and Pr denote the Reynolds number and 
Prandtl number, respectively. Then the governing 
equations may be written 2 


dn m dfm 

dt dz 


d 2 f„ 


or 


d z 2 
V-/L = 0 


= 0, m = 1,2,3, (4.10) 
(4.11) 


where 


v d M (A L\ 

v \dx’ at 


dz' dt) 


(4.12) 


and 


m=i - 2 ' 3 ' < 4i3) 


The integral form of (4.11) is then 

h m ds = 0, 


/ - 

JS(V) 


(4.14) 


where V is an arbitrary region in the space- time E%. 

Let Ei be discretized by solution elements 
SE(j, n) and conservation elements CE(j,n) of Type 
I, as described in [2] and [12]. Then, associated with 
each SE(j,n), let u m , / m and f m be represented in 
discrete form through the linear expansions 

t t* m (z,t\j,n) W (4.15) 

(«m)J + (*“*") 

fn(z,t-,j,n) d = (4.16) 

(/m)" + (* *i) + {fmt) j {t — t n ) 

n) d u (4.17) 


(fm)j + ( fmx)j ( x ~ x j) + ( fmt ) j (t O 

m = 1,2,3. 

The discrete coefficients (ii m )”, (tw)?, and 
m = 1>2,3, represent nine unknowns associated 

with each SE(j, n). The coefficients (/ m )", (/m)"> 
etc., are known functions of the (««»)". (u»i*)"> an< ^ 

(«mi)"- By requiring h' m d = (/A “ ^f-> u m)> m = 
1,2,3, to satisfy (4.11), each («mt)" may also be ex- 
pressed as a known function of the (um)" and (« mr )" . 
Thus, the marching problem from one time level to the 
next may be expressed entirely in terms of the discrete 
unknowns (w m )" and («mr)"- The marching condi- 
tions for and (u mr )" are obtained by requiring 

h' m to satisfy 


i 


S{CE(j,n) 




ds = 0, 


(4.18) 


for all (i,n). The reader is referred to Reference [2] 
for the details. 

In Figure 23, numerical solutions are presented to 
the shock tube problem with viscosity at a Reynolds 
number of 12,000. Note that the shock is resolved in a 
single mesh point, and that the contact discontinuity 
is spread over only three mesh points. No flux limiters, 
weighting functions, or other ad hoc parameters were 
used in the calculation. 

Summary 

A new numerical method is being developed for 
solving the Navier-Stokes equations. In this paper, the 
authors’ 2-D, steady-state implicit scheme has been 
generalized and extended to nonorthogonal Cartesian 
grids. Calculations of a variety of flow fields have 
demonstrated the ability of the scheme to accurately 
resolve developing boundary layers on coarse grids. In 
addition to its high accuracy, the main advantages of 
the present scheme are its conceptual simplicity, rig- 
orous enforcement of both the integral and differential 
forms of the governing equations, and its ability to 
be extended to three dimensions in a straightforward 
manner. Work in progress is directed toward extend- 
ing the scheme to three dimensions and developing a 
time-iterative solution technique for solving the dis- 
crete system of nonlinear equations. Finally, a scheme 
which uses cubic expansions to respresent the discrete 
primitive variables is also under development. 


11 



References 

[1] Scott, J.R. and Chang, S.C., “A New Flux Con- 
serving Newton’s Method Scheme for the Two- 
Dimensional, Steady Navier-Stokes Equations,” 
NASA TM 106160, June, 1993 (Accepted for pub- 
lication in the International Journal of Computa- 
tional Fluid Dynamics). 

[2] Chang, S.C., “The Method of Space-Time Con- 
servation Element and Solution Element - A New 
Approach for Solving the Navier-Stokes and Euler 
Equations,” Submitted to the Journal of Compu- 
tational Physics . (First published as “New De- 
velopments in the Method of Space-Time Conser- 
vation Element and Solution Element - Applica- 
tions to the Euler and Navier-Stokes Equations,” 
NASA TM 106226, August, 1993.) 

[3] Scott, J.R., “A New Flux-Conserving Numerical 
Scheme for the Steady, Incompressible Navier- 
Stokes Equations,” NASA TM 106520, April, 
1994 (Submitted to the International Journal of 
Computational Fluid Dynamics). 

[4] Chang, S.C. and To, W.M., “A New Numerical 
Framework for Solving Conservation Laws — The 
Method of Space-Time Conservation Element and 
Solution Element,” NASA TM 104495, August, 
1991. 

[5] Anderson, D.A., Tannehill, J.C., and Pletcher, 
R.H., Computational Fluid Mechanics and Heat 
Transfer, (Hemisphere, 1984), pp. 190-193. 

[6] Srinivasan, K. and Rubin, S.G., “Segmented 
Multigrid Domain Decomposition Procedure For 
Incompressible Viscous Flows”, International 
Journal For Numerical Methods in Fluids , Vol. 
15, pp. 1333-1355, 1992. 

[7] Schlichting, H., Boundary-Layer Theory , 7th. ed., 
McGraw-Hill Series in Mechanical Engineering, 
Holman, J.P., Consulting Ed., (McGraw-Hill, 
1979), pp. 135-140. 

[8] Dill, L.H. and Himansu, A., “Development and 
Validation of Algorithms for the Incompressible 
Navier-Stokes Equations Based upon the Space- 
Time Solution (STS) Element Method”, NASA 
Contractor Report NAS3-27014. 

[9] Messiter, A.F., “Boundary-Layer Flow Near The 
Trailing Edge of a Flat Plate”, SIAM J. Appl. 
Math., 18, 1, 241 - 257. 


[10] Sfcewartson, K., “On the flow near the trailing 
edge of a flat plate”, II, Mathematika, 16, 106 
- 121 . 

[11] Scott, J.R., Dill, L.H., Oyediran, A., and 
Greenspan, B., “A New High-Order Numerical 
Method for Solving Ordinary Differential Equa- 
tions”, NASA TM in preparation. 

[12] Chang, S.C. and To, W.M., “A Brief Descrip- 
tion of a New Numerical Framework for Solving 
Conservation Laws — The Method of Space-Time 
Conservation Element and Solution Element,” 
Proceedings of the Thirteenth International Con- 
ference on Numerical Methods in Fluid Dynamics , 
Rome, Italy, 1992, Napolitano, M. and Sabetta, 
F., eds., Lecture Notes in Physics 414, Springer- 
Verlag. Also published as NASA TM 105757. 


Appendix 

Local Constraints; 

Conservation of Mass 


**i,o + V 0>1 = 0 

(A.1) 

2«a,o + Vi.i = 0 

(A.2) 

2 t> 0>a + **i,i = 0 

(A.3) 

Conservation of x-Momentum 


**1,0 **0,0 "h **0,1 *\>,0 Pl,0 

(A.4) 

- ^[(2 **0.0 + *-,,) + |(4« 3 ,o-V,,i)] = 

0 

2 (2 t* 0 ,o **i,o + **?,„ + p 3 ,o) + ** 1,1 Vo, o 

(A.5) 

+ *>1,1 **0.0 + **1,0 *>0,1 + **0,1 Vl,0 = 0 


2(**o.o«o.a + V 0 ,o **o,3 + fo,i**o,i) 

(A.6) 

+ 2 t* Xil ti 0 .o + 2« l , 0 «o,i + Pi,i = 0 


Conservation of y-Momentum 


Uo ,0 V l ,0 d" ^ 0,0 Po,l 

(A.7) 

- •^[(2v 3 . # + *1i.i) + ^ (4 t>„, 3 = 

0 
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(w 0>0 ty a W 0t o ^ 0 , 3 ) 

(A.19) 

2 (tto ,0 ^ 3,0 4" ^ 0,0 ^3,0 4 ti^o) 

(A. 8 ) 

*” (Wo,o W a ,o W 0| o ^ 3 , 0 ) — 0 


4 2d 1i1 v 0i0 + 2 v 1>0 v 0tl + Pi.i = 

0 

Additional Eliminated Variables: 




Pa , 0 = Wo,o Wq,3 W 00 W a>0 

(A. 20) 

2 (2 v 0i0 ty 3 + *4 Po.a) + w 1 ( 1 1> 0|0 

(A.9) 

~ 2 ^^*° U °' 1 Vl ’ 0 ^ 


-4 v lfl ti 0>0 -4 w 1 > 0 1? 0|1 + w 0fl v lt0 = 

: 0 



General Eliminated Variables: 


Pl,l = 2 (t4o,0 W 0 ,2 W o ,0 Wo, 3 ) 

(A.21) 

Vo.i = “ W 1(0 

(A.10) 

Flux Vector Comnonents: 


v lt i — — 2 u 3 ,o 

(A.ll) 

Conservation of Mass 


ti t ,i — 2 v 0 , a 

(A.12) 

A*. = Wo.o 

(A .22) 

Pl.O = 

(A.13) 

Af l0 = t*i,o 

(A.23) 

2 

— til .0 tio ,0 W 0 ,i t> 0,0 ”4 ^ (Wa ,0 "4 

Wo, 2 ) 

a:., = Wo.i 

(A .24) 

Po,l = 

(A.14) 



2 

-Vi.otio.o -4 1 i li 0 ty<> + ^(^.0 + 

v 0 ,j) 

A*o = w».o 

(A.25) 

Po ,3 = 

(A.15) 

A*, = -2 hi, = — 2v 0>a 

(A .26) 

1/2 . 

Wo ,0 W a> o Wo ( o Wq,3 2 (« 1,0 ”4 Wq ,1 

Wi.o) 

A*2 = Wo, 2 

(A.27) 

Constant-Pressure- Gradient Scheme: 



(A .28) 

Local Constraints: 


Afl.o = w 0 , 0 

Pa.o = Vq,o Wo, 3 W 0f o W a ,o 

(A.16) 

A?,o = Wx.o 

(A.29) 

1 

'cT 

M tO 
e 

+ 

c 

0 

M 

c; 

0 

II 

0 




L 


A?., = — A* 0 = — «i,o 

(A .30) 

P\ t x = 2 (w 00 U 0f3 W 0>0 W 0j2 ) — ^ 

(A.17) 

A» 0 = v,.o 

(A.31) 



A?,2 = -2Af i0 = -2u 3 . 0 

(A.32) 

|>i,i = 2 (ty 0 W 2 ,o W 0>0 V 3 ,o) = 0 

(A.18) 

a :, 2 = wo , 2 

(A.33) 


General Scheme: 

Local Constraint: Conservation of z- Momentum 
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2 

2 

(A.34) 

*'2 = «0.oW 0 .3 + ^ 0,0 1 * 0,2 - «l,O« 0 ,l 

(A. 45) 

<« + Po.o - 





Conservation of y-Momentum 


*f,. = -*k 

= 

(A .35) 

* 0,0 = u 0.0 t'o.o - ^(« 0.1 + v l>0 ) 


+ «0.0«l,0 + 

2 , 

Jfe ( “ oa 

“ w a,o) 

(A.46) 


A* x — 2 ti 0t0 ti 0> i 4 v 0 ,o w i,o tio.o Vi.o (A.36) 


4 ® v o,a) 


Ko = «o.c v i,o + V 0 ,o «i.o “ ^(*> 2,0 " t>o,a) (A.47) 


^ 0,1 v OfO u o,i u o,o w i,o Re^ U ° 2 (A.48) 


h* Q = = (A.37) h* 20 = Uo.oi'j,# + v 0t0 « 3(0 4- v l>0 « lf0 (A.49) 


Wo^o u a,o 4* Vo .0 v o,s 4 2 u *.° 2 t4 °* 1 Vl *° 


K tl = — 2 h* 2 = 


(A.38) 


— 2 (ti 0 ,0 V 0 ,3 4 Vo,0 U 0,3 U l,0 U 0,x) 

h * t 2 ~ 2 ti 0(0 t< 0|3 + ti 0t0 ti 3>0 v 0fQ v o,a (A.39) 

2 1 2 1 
4 ti 0>l - -u l 0 - -U 0tl v lf0 

Aj.o = «0, 0*0,0 - ^(«o,i 4 V 1>0 ) (A.40) 


= ti 0 , 0^1,0 4 t>0,0«1.0 - *^(*3,0 ” *0,2 ) (A.41) 


i — *o,o **o,x u 0|0 ti lf0 ReS^*' 2 (A.42) 


h y 0 = « 0f0 v 3f0 4 v 0t0 u 2t0 -I- t; 1(0 « lt0 (A.43) 
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Figure 4. Channel geometry and discretization. 



Figure 3. Cell orientation for flux expressions. 


Figure 6. Developing channel flow at 1, 3, 5, and 11 
channel heights downstream. Re = 100. 39 x 6 grid. 
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Figure 10. Comparison of v velocity component with 
Blasius solution. Re = 100,000. 110 x 28 grid. 



Figure 11. Pressure coefficient in the trailing edge 
region. Re = 100,000. 110 x 22 grid. 







Figure 13. Comparison of numerical results from the x 

generalized scheme and the constant-pressure-gradient 

scheme. Channel flow at Re = 2000. z = 1> 5, 9. Figure 16. Nonorthogonal Cartesian grid. 



Figure 14. Comparison of numerical results from the 

generalized scheme and the constant-pressure- gradient Figure 17. Local coordinate system for flux evalua- 
scheme. Pressure coefficient in the trailing edge region tion. 
of a flat plate airfoil. Re = 7500. 



Figure 15. Comparison of numerical results from the Figure 18 a. Computational grid for a converging 

generalized scheme and the constant-pressure-gradient channel. Ramp angle = 10 degrees. 32 x 8 grid, 

scheme. Pressure coefficient in the trailing edge region 
of aflat plate airfoil. Re = 1000. 
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Figure 18 b. Developing flow in a channel with a 
converging section. Re = 100. 


Figure 20 a. Computational grid for a diverging chan 
nel. Ramp angle = 5 degrees. 32 x 8 grid. 


Figure 19 a. Computational grid for a channel with a 
15 degree ramp. 32 x 10 grid. 
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Figure 20 b. Developing flow in a diverging channel 
with a 5 degree ramp. Re = 100. 



Figure 19 b. Developing flow in a channel with a 15 
degree ramp. Re — 200. 


Figure 21. Comparison between tt 0(1 and the first y 
derivative of the Blasius solution, z = .805. 110 x 28 
grid. Re = 100,000. 
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